home *** CD-ROM | disk | FTP | other *** search
/ FishMarket 1.0 / FishMarket v1.0.iso / fishies / 376-400 / disk_386 / xlispstat / src2.lzh / XLisp-Stat / ppnd.c < prev    next >
C/C++ Source or Header  |  1990-10-04  |  2KB  |  72 lines
  1. /* XLISP-STAT 2.1 Copyright (c) 1990, by Luke Tierney                  */
  2. /* Additions to Xlisp 2.1, Copyright (c) 1989 by David Michael Betz    */
  3. /* You may give out copies of this software; for conditions see the    */
  4. /* file COPYING included with this distribution.                       */
  5.  
  6. #include "xlisp.h"
  7. #include "osdef.h"
  8. #ifdef ANSI
  9. #include "xlsproto.h"
  10. #else
  11. #include "xlsfun.h"
  12. #endif ANSI
  13.  
  14. #define zero   0.0 
  15. #define half   0.5 
  16. #define one    1.0 
  17. #define split  0.42e0 
  18. #define a0     2.50662823884e0 
  19. #define a1   -18.61500062529e0 
  20. #define a2    41.39119773534e0 
  21. #define a3   -25.44106049637e0 
  22. #define b1    -8.47351093090e0 
  23. #define b2    23.08336743743e0 
  24. #define b3   -21.06224101826e0 
  25. #define b4     3.13082909833e0 
  26.  
  27. #define c0    -2.78718931138e0 
  28. #define c1    -2.29796479134e0 
  29. #define c2     4.85014127135e0 
  30. #define c3     2.32121276858e0 
  31. #define d1     3.54388924762e0 
  32. #define d2     1.63706781897e0 
  33.  
  34. /*
  35. c
  36. c   Algorithm as 111 Applied statistics (1977), vol 26 no 1 page 121
  37. c   Produces normal deviate corresponding to lower tail area of p
  38. c   the hash sums are the sums of the moduli of the coefficients
  39. c   they nave no inherent meanings but are incuded for use in
  40. c   checking transcriptions.  Functions abs,alog and sqrt are used.
  41. c
  42.  
  43. derived from AS111 fortran version
  44. */
  45.  
  46. double ppnd(p, ifault)
  47.     double p;
  48.     int *ifault;
  49. {
  50.   double q,r,ppn;
  51.  
  52.   *ifault = 0;
  53.   q = p - half;
  54.   if( fabs(q) <= split) {
  55.     r = q*q;
  56.     ppn = q *  (((a3 * r + a2) * r + a1) * r + a0) 
  57.             / ((((b4 * r + b3) * r + b2) * r + b1) * r + one);
  58.   }
  59.   else {
  60.     r = p;
  61.     if(q > zero) r = one - p;
  62.     if(r <= zero) {
  63.       *ifault = 1;
  64.       return(0.0);
  65.     }
  66.     r = sqrt(-log(r));
  67.     ppn = (((c3*r+c2)*r + c1) * r + c0) / ((d2 * r + d1) * r + one);
  68.     if( q < zero) ppn = -ppn;
  69.   }
  70.   return(ppn);
  71. }
  72.